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An understanding of the hydrophobicity of complex heterogenous molecular assemblies is crucial to charac¬ 
terize and predict interactions between biomolecules. As such, uncovering the subtleties of assembly processes 
hinges on an accurate classification of the relevant interfaces involved, and much effort has been spent on 
developing so-called “hydrophobicity maps.” In this work, we introduce a novel electrostatics-based mapping 
of aqueous interfaces that focuses on the collective, long-wavelength electrostatic response of water to the 
presence of nearby surfaces. In addition to distinguishing between hydrophobic and hydrophilic regions of 
heterogenous surfaces, this electrostatic mapping can also differentiate between hydrophilic regions that po¬ 
larize nearby waters in opposing directions. We therefore expect this approach to find use in predicting the 
location of possible water-mediated hydrophilic interactions, in addition to the more commonly emphasized 
hydrophobic interactions that can also be of significant importance. 
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I. INTRODUCTION 

The formation of complex mesoscopic structures is gen¬ 
erally driven by the self assembly of nanoscale build¬ 
ing blocks. Oftentimes, the assembly of individual units 
into a larger complex is driven by strong, direct inter¬ 
actions between the chemical constituents®*. In many 
cases, however, the assembly of molecular building blocks 
is triggered by the solvent and not by direct interac¬ 
tions between the individual units. For example, water- 
mediated interactions between nonpolar entities, known 
as hydrophobic effect ^El, induce the folding of proteins 
by minimizing the contact area between hydrophobic 
amino acid residues and water at the surface of the pro- 
tehPS. Both hydrophobic and electrostatic interactions 
are also thought to play a role in protein-protein binding 
processes 2 11 12 . Driven by such applications, many re¬ 
searchers have proposed empirical hydrophobicity scales 
for characterizing and predicting the relative hydropho¬ 
bicity of complex surfaces. In general, these approaches 
fall into one of two classes, surface-based and water-based 
hydrophobicity scales. 

Surface-based hydrophobicity scaleJ^®^ seek to pre¬ 
dict the relative hydrophobicity of a water-surface inter¬ 
face from the properties of the surface alone. Such map¬ 
pings have generally been used to predict the properties 
of proteins and other complex biological surfaces. Clas¬ 
sic hydrophobicity scales like the Kyte-Doolittle scale 13 
assign a value of relative hydrophobicity to each macro- 
molecular building block or residue. 

Another residue-based scale was introduced by Berne 
and coworkers to characterize the binding regions of pro¬ 
tein surfaces involved in protein-protein interactions, and 
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hundreds of surfaces were examined to predict if a dewet¬ 
tin g tr ansition could occur upon protein complex forma¬ 
tion^. This analysis indeed successfully uncovered sev¬ 
eral binding processes in which a dewetting transition 
was likely to occur, as confirmed by molecular dynamics 
simulations with explicit water molecule^*. 

More recent approaches have demonstrated that in¬ 
corporation of atomic-level details may provide an even 
better characterization of surfaces than classifying an en¬ 
tire residue as hydrophobic or hydrophilic. In particu¬ 
lar, Kapcha and Rossky (KR) have exploited the assign¬ 
ment of partial charges in atomistic models of proteins 
to classify individual atomic units as hydrophobic or hy¬ 
drophilic on a binary scald 16 . Predictions from these KR 
maps are in much better agreement with more detailed 
results than those from residue-based methods and have 
been used to characterize a wide range of protein surfaces 
and binding pocketd^. 

Despite these predictive successes, surface-based hy¬ 
drophobicity scales only implicitly account for the influ¬ 
ence of surface topography and the chemistry of nearby 
atomic units on the hydrogen bond network of water. 
In light of these limitations, new methods have been 
developed for providing a direct characterization of the 
nonlocal response of water to complex material surfaces. 
In particular, recent work has focused on the nature of 
density fluctuations in bulk water and at aqueous inter- 
facepES. 

Density fluctuations are intimately related to solvation 
thermodynamics through potential distribution theory, 
which can be used to relate the solvation free energy of 
a hard object to the probability of observing a cavity 
or volume of equal shape and size being empty 24 27 . A 
large cavity or extended hydrophobic surface disrupts the 
hydrogen bond network and density fluctuations are en¬ 
hanced relative to those in the bulk or near a hydrophilic 
surface. This disruption of the network is the key physi- 





2 


cal feature characterizing large scale hydrophobic effects, 
and this makes it easier to create a cavity near a hy¬ 
drophobic surface than a hydrophilic one. Previous work 
has exploited this fact, and used the solvation free en¬ 
ergy of hard shaped objects or cavities near a surface as 
a measure of its hydrophobicit}) 17 23 . 

Because this approach focuses on the response of water 
to a substrate, density fluctuation-based mappings can 
account for both chemical and topographical complexi¬ 
ties of a surface. However performing such calculations, 
especially for large asymmetric volumes, can be computa¬ 
tionally intense, often requiring several simulations with 
advanced umbrella-sampling techniques to obtain the sol¬ 
vation free energy of a single volume. Moreover, this type 
of mapping depends on the size, shape, and position of 
the probe volume to be emptied and in the absence of 
a general theory, different choices can lead to ambigu¬ 
ities in the estimate of hydrophobicity. While fluctua¬ 
tion methods very successfully discriminate between hy¬ 
drophobic and hydrophilic surfaces, little distinction is 
seen between surfaces with similar hydrophilicity but dif¬ 
ferent chemistries. 

In this work, we introduce a simple approach to char¬ 
acterize chemically and topographically complex materi¬ 
als based on the collective, long-wavelength electrostatic 
response of water to such surfaces. Excluded volume con¬ 
straints and the partial surface charges exploited in the 
KR scale can produce distortions and disruptions of the 
hydrogen bond network. Extended hydrophobic regions 
can break hydrogen bonds, while polar regions can lo¬ 
cally distort and pin bonds in the network. This compli¬ 
cated nonlocal response rearranges the molecular dipoles 
in water and is encoded in the electrostatic potential. 

As we will see, the long-wavelength component of 
the induced potential readily distinguishes between hy¬ 
drophobic and hydrophilic surfaces, and also provides ad¬ 
ditional insight into the net polarization of water at a 
polar interface. This is used to differentiate between hy¬ 
drophilic surfaces with similar contact angles but differ¬ 
ent chemical structures. We anticipate that this compu¬ 
tationally simple method of analysis, which requires only 
a single equilibrium simulation, will find use as a tool to 
efficiently characterize complex surfaces with nanoscale 
chemical and topographical patterning, like the surface 
of proteins involved in complex formation. However, the 
forces generating the distortions of the hydrogen bond 
network are by no means purely electrostatic, and direct 
physical connections to macroscopic measurements of hy¬ 
drophobicity like the contact angle or the solvation free 
energy of cavities will require more general treatments. 


II. ELECTROSTATIC RESPONSE AT UNIFORM 
PLANAR SURFACES 

We first introduce our approach and its underlying 
concepts through the study of uniform planar surfaces 
with varying hydrophobicity. In particular, we study the 
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FIG. 1. (a) Bare charge densities provide little information 
about the collective electrostatic response of water to model 
silica surfaces of varying polarity, £. (b) Upon Gaussian- 

smoothing, however, a long-ranged, collective dipole layer 
at the interface is apparent. This dipole layer character¬ 
izes the response of water to the polarity of the surface, (c) 
The Gaussian-smoothed electrostatic potential due to water 
provides a slowly varying metric quantifying this interfacial 
dipole layer, (d) The difference A Vs in the long-ranged elec¬ 
trostatic potential in bulk water and in the effective vacuum 
inside the silica surface, or at the surface itself as indicated by 
the dashed line, characterizes the nonlinear response of water 
to the magnitude of the dipole of the surface. 


response of water to the model atomically-detailed and 
corrugated silica-like surfaces of Giovambattista, Rossky, 
and Debenedetti 28 29 . These artificial surfaces are dipo¬ 
lar, with the dipole moment of surface groups created by 
negative partial charges placed on the uppermost layer of 
oxygen sites and positive charges on the subsequent layer 
of silicon atoms. In order to tune the relative polarity of 
the surface, the magnitude of the atomic charges, g, is lin¬ 
early coupled to a parameter £, such that the dipole mo¬ 
ment of each surface unit is given by £p = (0, 0, £gdosi)> 
where dosi is the oxygen-silicon bond length. A nonpo¬ 
lar, hydrophobic surface with a contact angle of 6% « 109° 
is obtained when £ = 0. Increasing the magnitude of £ 
leads to an increase in the hydrophilicity of the surface, 
as characterized by a decrease in the contact angle 
and there is a slight asymmetry with respec t to th e sign 
of £. For example, 0\ ~ 90°, while 0_i « 85 ^ 28 * 29 l 

In Figure [l^i, we show p q (z ), the variation in the z- 
direction of the bare charge density of water, 
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qi and position r$(R) in configuration R, and (• • •) in¬ 
dicates an ensemble average over all configurations and 
we have averaged over the remaining two coordinates in 
r. Because of the rapidly varying molecular scale sur¬ 
face structure and charge distributions, little insight into 
the collective response of water to the surface is garnered 
from p q (z ), and an alternative description of the ordering 
of inter facial water is needed. 

A clearer measure of the long-wavelength electrostatic 
response of water induced by a heterogeneity is given by 
the Gaussian-smoothed charge density p qcF (?) that natu¬ 
rally arises within the context of the local molecular field 
(LMF) theory of electrostatic^!. Here 

P qr7 ( r) = J dr'p g (r')p G (\r - r'|), (2) 

Pg( r) is a normalized Gaussian function given by 


Pg( r) 



(3) 


and a is a smoothing length used to average over 
short-ranged rapidly varying components of the electro¬ 
static interactions. Previous work has shown that use¬ 
ful smoothing results wheji a is chosen on the order of 
nearest-neighbor distance^H^, and we choose a = 4.5 A 
herein. This type of smoothing over molecul ar siz es, often 
discussed in elementary electrostatic text d 33 * 3 ^ !, washes 
out the effects of atomic-scale details, and it has been 
argued that p qa ( r) provides a good qualitative charac¬ 
terization of the underlying long wavelength electrostatic 
response of fluids 31 32 . Note that in the present case when 
the bare charge density p q ( r) is known directly from sim¬ 
ulation, LMF theory is used only to motivate the con¬ 
sideration and usefulness of the smoothed charge density 
p qa ( r). 

Indeed, as shown in Figure [TJ d, the Gaussian-smoothed 
charge densities readily indicate the formation of an in¬ 
terfacial dipole layer arising from subtle changes in the 
orientation of water molecules near the surface as they 
are increasingly pinned by the partial charges of the 
model surface. The formation of a dipole layer is appar¬ 
ent even at the nonpolar surface with £ = 0, as detailed 
in previous wor M—lSMgej 

As £ is increased from zero in the negative direction, 
this interfacial dipole grows in magnitude. The dominant 
molecular orientation producing the interfacial dipole for 
£ <C 0 corresponds to that sketched in Figure i=. Here 
a non-negligible fraction of interfacial water molecules 
point one O-H bond directly toward the surface, and 
form hydrogen bonds with the negatively charged “sil¬ 
icon” atoms. The collective interfacial dipole at £ = 0 
arises from the smaller fraction of water already exist¬ 
ing in this orientation, although most water molecules at 
the interface then have an orientation where their H-O- 
H plane is parallel to the surface itself, which makes no 
contribution to the dipole layei!^. 


Conversely, as £ is increased from zero in the posi¬ 
tive direction, the magnitude of the interfacial dipole de¬ 
creases and eventually changes sign. This is consistent 
with interfacial water molecules pointing an O-H bond 
away from the interface with increasing probability as £ 
is increased, as evidenced by the probability distributions 
P(0 O u) of the angle made by the OH bond vector of a wa¬ 
ter molecule and the surface normal shown in Figure^. 

Although the Gaussian-smoothed charge density has 
provided a physically suggestive description of the re¬ 
sponse of water to a uniform planar surface, a more quan¬ 
titative metric is needed to characterize patterned and 
corrugated surfaces and protein complexes. To that end, 
we focus attention on Vs(r), the slowly-varying, long- 
wavelength portion of the electrostatic potential due to 
the solvent. This is related to the smoothed charge den¬ 
sity p qa ( r) by Poisson’s equation 

V 2 V s (r) = -p^( r), (4) 

with notation chosen to be consistent with previous 
work 30 . This potential, shown in Figure [l}l , smoothly 
transitions from zero in the vacuum region well inside 
the silica surface to a value of A Vs in the bulk region. 
Note that trends in A Vs can be evaluated using other 
reasonable limits, like the difference between the bulk po¬ 
tential and the dashed line in Figure [1J at the silica sur¬ 
face, without qualitatively affecting the characterization 
of the interface. We use this fact below when generating 
an electrostatics-based map of a water-protein interface. 

We propose that AVs can be utilized as a quantitative 
measure of the response of water to complex surfaces. In¬ 
deed, AVs(£) shown in Figure [l]i quantifies the nonlinear 
response of the mesoscopic surface dipole layer to changes 
in £, such that the magnitude and sign of water polariza¬ 
tion at the surface is captured. Additionally, the relative 
differences in AVs(£) are due exclusively to structural 
rearrangements near the surface (dipolar contributions), 
because the constant quadrupolar (Bethe potential) con¬ 
tribution to the potential is the same for all £-valued^. 

Moreover, recent work has emphasized that polar hy¬ 
drophilic surfaces, which locally pin water molecules, are 
fundamentally different from nonpolar hydrophilic sur¬ 
faces that have artificially large Lennard-Jones-like at¬ 
tractions. These latter surfaces simply pull interfacial 
waters closer to the surface, increasing the contact angle, 
but do not introduce significant chan ges in water struc¬ 
ture about the position of this interface 39 ! Therefore, the 
local pinning of interfacial waters through polar interac¬ 
tions like hydrogen bonding underlies the hydrophilicity 
of realistic surfaces, and the electrostatic mapping intro¬ 
duced here captures this behavior. 

Finally, we note that in the simple case of a system 
with a slab-like geometry, the potential difference ob¬ 
tained from the one-dimensional Vs (z), AVs, is exactly 
equal to that estimated from the bare electrostatic po¬ 
tential V(z), AV. This equality is due to the fact that 
Gaussian smoothing of the potential conserves the first 
two nonzero multipole moments of a charge distribution 
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FIG. 2. (a) The probability distribution of the angle made by the surface normal and the OH bond vector, P(0oh), illustrates 

that water polarizes differently in response to surfaces with positive and negative £. Gray arrows indicate the direction of 
increasing £. (b) For £ ^ 0, interfacial water molecules tend to point an OH bond toward the bulk, while (c) interfacial waters 
point an OH group toward the surface for $ « 0. A small portion of each surface of varying £ is indicated by the colored 
rectangle, while the oxygen site is colored red and the hydrogens are colored light gray. 



FIG. 3. (a) An otherwise apolar silica surface, composed of oxygen (green) and silicon (yellow) atoms, is patterned with 

a single, hydrophilic hydroxyl group, shown in red (oxygen) and white (hydrogen), (b) The potential difference AVs(x,y) 
(in Volts) from the three-dimensionally Gaussian smoothed charge density reveals that the polar and nonpolar portions of the 
surface polarize nearby water differently, enabling the mapping of surfaces with complex chemistries, (c) The characterization of 
this patterned surface provided by the solvation free energy of a hard cuboidal solute, /3A/jl(x, y) (see text for details), illustrates 
that the electrostatics-based map is consistent with previous, water density fluctuation-based surface mapping techniques. A 
green diamond indicates the position of the hydroxyl group oxygen in (b) and (c). 


(see the Appendix), the dipole and quadrupole in the case 
of water, and the potential difference across two phases in 
a slab geometry depends only on these two moments for 
nonionic systems^®. However, this is true only in the 
special case of slab-like symmetry. Gaussian-smoothing 
over molecular length scales is critical for quantifying the 
electrostatic response to complex surfaces, where Vs(r) 
is three-dimensional and has no apparent symmetry, as 
we discuss below. 


III. CHARACTERIZATION OF CHEMICALLY 
PATTERNED SURFACES 

The above analysis can be readily generalized to more 
complex, patterned surfaces with less symmetry. Here 
we show that the potential difference as a function of 


lateral position, AVs(x,y), can be used to characterize 
the long-wavelength perturbations of the H-bond network 
induced by local changes in the chemical patterning of 
a planar molecular surface. Gaussian smoothing of the 
three-dimensional charge density in Eq. © ensures that 
relevant non-local perturbations of the water structure 
due to patterning of the surface in the xy -plane are cap¬ 
tured, in addition to those in the ^-direction, parallel to 
the surface normal 

We first consider a purely hydrophobic (uncharged) sil¬ 
ica surface, and add a single hydrophi lic sit e, consisting 
of a hydrogen-bonding hydroxyl groupF^ES as shown in 
Figure |3 ^l. A clear picture of the electrostatic response 
of water emerges upon examination of A Vs (x,y). Long- 
wavelength perturbations of the water H-bond network 
extending over roughly 10 A in the plane of the surface 
are found when a single hydrophilic site is added to the 
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otherwise hydrophobic surface, Figure ^jp. This is consis¬ 
tent with the concept that a single hydrophilic site can 
substantially pin water in its vicinity 21 4 -^, and also is in 
accord with the results detailed above for homogeneous 
surfaces. As the number of hydrophilic sites is increased 
to three and seven sites, an increasingly larger region of 
water is perturbed, as illustrated in Figures |4^i and [4 |d, 
respectively. 

We also consider surfaces complementary to those in 
Figures [3] and [4j created by replacing hydrophilic groups 
with hydrophobic groups, and vice versa. For example, 
the complement of the surface in Figure [3^i has a sin¬ 
gle apolar unit at the center of a surface covered in hy¬ 
droxyl groups. The electrostatics-based maps for these 
hydrophilic surfaces with one, three, and seven hydropho¬ 
bic sites are shown in Figure [5^t, [5)3, and [5J3, respec¬ 
tively. Our findings are consistent with the idea that 
a hydrophobic site does not perturb the H-bond network 
near a hydrophilic surface nearly as much as a hydrophilic 
site does on an otherwise hydrophobic surfac^ 21 . This 
can clearly be observed by comparing Figure [3)3 to Fig¬ 
ure UK where a single hydrophobic site hardly perturbs 
the structure of water at the surface. 

We additionally note that the response of water to a 
small hydrophobic patch on a hydrophilic surface is not 
equal to that above a uniform hydrophobic surface (a 
similar statement holds for the complementary surfaces). 
This equality becomes true for large enough patches, but 
the collective effects due to maintaining the H-bond net¬ 
work across the patch boundaries restricts the orientation 
of inter facial water molecules, and leads to larger values 
of A Vs (x,y) above these small patches than is expected 
from uniform surfaces with the same chemistry, see Fig¬ 
ure [lji. 

We further compare our electrostatics-based mapping 
procedure with a density fluctuation-based mapping in 
Figure^. The fluctuation maps shown here were gener¬ 
ated by calculating the probability P V (N ; x, y) of observ¬ 
ing N water molecules in a probe volume v = a x a x 4 A 
centered at the first peak of the nonuniform density of 
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FIG. 4. Water molecules at the surface are significantly 
polarized when (a) three and (b) seven hydrophilic (hydroxyl) 
groups are added to an otherwise hydrophobic model silica 
surface, as indicated by the AVs(x,y) plots shown here (in 
Volts). The positions of the hydroxyl group oxygen sites are 
indicated by the green diamonds. 


water at the surface for all x and y. The free energy of 
emptying this volume over the surface is then given by 
(3A/j,(x, y) = — In P v (0; x, y) and this is used to probe the 
relative hydrophobicity of the surface. 

This fluctuation-based mapping provides a description 
of interfacial water in accord with the electrostatics-based 
approach presented herein. For all surfaces studied here, 
consistency is found between the fluctuation-based and 
electrostatic-based maps of relative hydrophobicity. 

Although both fluctuation-based and electrostatics- 
based maps discriminate well between hydrophobic and 
hydrophilic regions of a surface, the electrostatics-based 
mapping provides additional useful information about 
the polarization of water at the hydrophilic portions of 
the surface that standard fluctuations-based methods do 
not resolve. To illustrate this point, we construct a model 
surface with regions of differing polarity shown in Fig¬ 
ure [6^1. The first patch (red) has four surface units with 
£ = 1, while the other two (blue) correspond to units 
with £ = — 1, one with three units and one with two. We 
expect that water in the vicinity of the surface will be po¬ 
larized in opposite directions above the different dipolar 
units. The remainder of the surface is nonpolar (gray). 

We find that AVs(x, y) can readily distinguish between 
the regions of differing polarity. Water in the vicinity of 
the £ = 1 patch responds in a manner that yields a posi¬ 
tive A Vs, while a negative A Vs is obtained over £ = — 1 
regions. Surprisingly, the signs of these potentials are op¬ 
posite to those obtained at uniform surfaces, stemming 
from the small size of the patches and the fact that most 
of the surface is apolar. Water molecules above the apo¬ 
lar surface do not penetrate into the grooves of the sur¬ 
face, as do waters near uniform £ = 1 and £ = — 1 sur¬ 
faces. Penetration of water into the surface grooves near 
the small patches in Figure |6^i would lead to significant 
distortions of the hydrogen bond network, which are not 
compensated by these weakly hydrophilic patches. In¬ 
stead, the orientation of water molecules in the vicinity 
of these patches is dominated by the tendency for their 
dipoles to align with the dipolar field of these patches. 

This effect is due to the weakly hydrophilic nature of 
these artificial patches and the collective nature of the 
water hydrogen bond network. We expect a transition in 
A Vs with patch size, until the infinite limit shown in Fig¬ 
ure EE is reached. Indeed, differences in AVs with patch 
size are already observed in moving from two to three 
dipolar units (top right to central blue patches in Fig¬ 
ure UK)’ the magnitude of AVs increases. Such context- 
dependent behavior cannot be captured by focusing on 
the properties of the surface alone. In this case, surface- 
based mapping would predict no difference between these 
patches and the uniform surface because they have the 
same basic charges, despite the very different response of 
water based on the environment. These results provide a 
dramatic illustration of both the flexibility and subtlety 
of the nonlocal response of the hydrogen bond network to 
various perturbations, and the ability of the electrostatic 
maps to resolve their effects. 
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FIG. 5. The electrostatics-based mapping provided by AVs(^,y) (V) additionally provides a characterization of the response 
of water to hydrophilic silica surfaces with (a) one, (b) three, and (c) seven hydrophobic, apolar units, the positions of which 
are indicated by the diamonds. 



FIG. 6. (a) Snapshot of the model surface where the dipole moment of the surface units in the polar regions colored red and 
blue have opposite sign, while the rest of the surface is apolar (gray), (b) Potential difference AVs(^,y) due to the Gaussian- 
smoothed charge density above the patterned silica surface, (c) AVs(^,y) obtained above the inverse of the surface in panel 
(a), where the signs of the dipoles on the polar regions have been switched. Note that the opposite surface dipole moments 
induce differing responses in nearby water molecules, as captured by the smoothed electrostatic potential difference shown here. 


We also consider the “inverse” of the surface in Fig¬ 
ure UK which has the same pattern, but the signs of 
the dipoles are opposite to those of the surface in Fig¬ 
ure it- Thus, we expect water in the vicinity of these 
patches to be polarized in the direction opposite to those 
of the inverse surface. Indeed, the electrostatic mapping 
described by A Vs {x,y) demonstrates this to be true, as 
shown in [6^. This illustrates that the electrostatic map¬ 
ping of interfaces readily uncovers the qualitatively dif¬ 
ferent polarization of water at polar surface units, in ad¬ 
dition to the relative hydrophobicity of regions of the 
surface. 


IV. MAPPING PROTEIN SURFACES 

The surface of the protein Hydrophobin II (HFBII) has 
both a hydrophilic and a hydrophobic region, allowing 
the protein to effectively act as a biomolecular Janus par¬ 
ticle and assemble at hydrophobic-aqueous interfaces 44 . 
HFBII has been analyzed using classic hydrophobicity 


scales, in addition to the water density fluctuation-based 
techniques described in the Introductio n 18 1 21 . The wide 
array of chemical and topographical complexities on the 
HFBII surface thus presents an ideal test case for ex¬ 
tending the electrostatic mapping technique to protein 
surfaces. 

In contrast to simple planar interfaces, one cannot gen¬ 
erally define a vacuum region on the protein side of the 
protein-water interface. Instead, we can define a water- 
protein interface, s, and evaluate the difference between 
the value of the long-ranged potential in the bulk and 
at this interface, AVs(s). This modified metric can still 
distinguish between regions of varying polarity, since our 
qualitative findings for planar interfaces would be unal¬ 
tered by shifting the vacuum reference to the position of 
the interface (z « 0 A in Figure [I] for example). 

In the results shown here, we follow previous work= 
and evaluate AVs(s) at the Willard-Chandler smoothed 
interface defined by the protein heavy atoms alond^. We 
choose a low value of the density to define the interface, 
so that this surface envelopes the heavy atoms of the 
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FIG. 7. (a,c) Two views of the Kyte-Doolittle residue-based 

mapping of Hydrophobin II (HFBII); increasing hydropho- 
bicity coincides with the darkness of the residue. Panel (a) 
highlights large hydrophobic patch on the surface of the pro¬ 
tein, while panel (c) focuses on hydrophilic regions. (b,d) 
The analogous two views of the electrostatic mapping of the 
water-HFBII interface. The protein surface as described by a 
Willard-Chandler smoothed interface for the protein heavy 
atoms is shown as a meshgrid and is colored by value of 
the AVs(s). Red represents negative values of the poten¬ 
tial, blue indicates positive values of the potential, and white 
corresponds to values near zero, such that slightly negative, 
pale red/pink indicates values of the A Vs observed above hy¬ 
drophobic surfaces. Figures were made using VMD^l 

protein and is similar to a typical solvent excluded surface 
used commonly throughout the literature. As illustrated 
in Reference mi because of the Gaussian smoothing in 
Vs(r), many molecular details in the protein substrate 
are unimportant. Thus the surface mapping procedure is 
rather insensitive to the precise definition of the protein- 
water interface, and many other reasonable choices could 
be made as well. 

The electrostatics-based mapping of HFBII is com¬ 
pared to the classic Kyte-Doolittle hydrophobicity scale 
in Figure [7| The Kyte-Doolittle scale is indicated by 
the color of the residues in Figure [7^t and [TJd, such that 
hydrophobicity is proportional to the darkness of the 
residue color. The value of AVs(s) in Figure [7)3 and[7]i 
is indicated by the color of the grid points on the mesh, 
which corresponds to the Willard-Chandler surface of the 
protein. The average interface is colored such that red 
represents large negative values of AVs(s), blue indicates 
large positive values of the potential, and white corre¬ 
sponds to values near zero. With this scaling, values 
of A Vs (s) indicative of hydrophobic surfaces are colored 
pink/pale red. 

The electrostatic mapping procedure readily uncovers 
the large hydrophobic patch of HFBII shown in Figure [7^i 
and [7)3. The value of AVs(s) in the region of this patch 
is consistent with that for a hydrophobic surface. How¬ 


ever, note that there is a highly hydrophilic group at the 
edge of this hydrophobic patch (bottom left, red region). 
The polarization caused by this patch does not dissipate 
instantly, and the size of the hydrophobic patch is pre¬ 
dicted to be smaller than that from surface-based scales. 
This agrees rather well with findings of other water-based 
characterizations of HBFlJiHlij. 

As was detailed in the previous section, polar surfaces 
can polarize water in opposing directions. The maps of 
HFBII indicate the formation of such patterns in a bio¬ 
logical context. The blue (positive AVs(s)) portion of the 
map in Figure [7)3 is dominated by the response of water 
molecules to the positively charged lysine residue in the 
region. Analogously, the red (negative AVs(s)) region on 
the right of Figure [7)3 describes the response of water to a 
negatively charged aspartic acid residue. The identifica¬ 
tion of such hydrophilic regions of differing polarity may 
be crucial to understanding the binding of biomolecules. 
For example, these oppositely charged residues come into 
close contact with those of neighboring HFBII molecules 
when forming multimeric assemblies^, and a character¬ 
ization of the polarization of water near polar regions of 
proteins may help predict their involvement in biomolec- 
ular interactions^. 


V. CONCLUSIONS 

In this work, we have introduced a novel electrostatics- 
based method for characterizing the long-wavelength, 
collective response of water to complex surfaces. This ap¬ 
proach effectively coarse-grains over molecular scale de¬ 
tails in order to uncover the underlying electrostatic be¬ 
havior of the system, enabling one to distinguish larger 
scale rearrangements of the hydrogen bond network of 
water at an interface and consequently differentiate hy¬ 
drophobic and hydrophilic surfaces. This approach is 
quite efficient computationally, requiring only a single 
equilibrium simulation of a complex surface in explicit 
water. 

In addition to distinguishing between hydrophobic and 
hydrophilic regions, the approach introduced here can 
further distinguish the relative polarization of water in 
response to hydrophilic groups. An understanding of the 
relative polarization of water at surfaces is important for 
a range of processes, including understanding and pre¬ 
dicting ion adsorption to biomolecular interfaces. It has 
been postulated that the sign of the surface potential 
can drive ions to adsorb to interfaces, at least in classical 
modeld^. As such, if adsorption were to be favorable at 
a patterned surface like that in Figure [6j for example, 
one might expect that ions of opposite sign would adsorb 
to regions of opposite dipole moments. 

Differentiating the relative polarity of water in interfa¬ 
cial regions may also be useful in the description of bind¬ 
ing processes in general. In particular, it seems plausible 
that polarization-ba sed h ydration repulsion could occur 
between two surface^®^. For example, consider placing 






in close proximity two surfaces that have maps like those 
in Figure [ 6 J 3 , such that the blue regions align. As the 
two surfaces are brought closer together, polarized inter¬ 
facial water may “interfere deconstructively” and effec¬ 
tively repel the two surfaces, because the hydrogen bond 
network is polarized in different directions everywhere 
but the midpoint between the inter faces, where the po¬ 
larization must vanish by symmetry^ 51 . Therefore, to 
achieve small interplate distances, confined water must 
be (unfavorably) depolarized. 

In contrast, one might expect that bringing together 
surfaces with regions of opposing polarity, like the red 
and blue regions in Figure [ 6 j would instead lead to hydro¬ 
gen bond networks that interfere constructively , because 
the inter-surface water is polarized in the same direction 
everywhere between these patches. In this case, unfa¬ 
vorable depolarization may not occur, minimizing this 
contribution to the effective repulsion between surfaces 
and possibly leading to an effective attraction between 
hydrophilic surfaces 51 52 . 

However, in more general cases of molecular assem¬ 
bly, the polarization of water molecules, and therefore 
AVs(r), depends on a multitude of parameters includ¬ 
ing the separation and relative orientation of the confin¬ 
ing surfaces, and more than bulk equilibrium simulations 
of individual surfaces are required for a quantitatively 
accurate description of such phenomena. Nevertheless, 
the electrostatic mapping described herein should pro¬ 
vide qualitative insight into many relevant features of 
solvent polarization effects needed to understand hydra¬ 
tion mediated forces. 


VI. METHODS 
A. Simulation of Silica Surfaces 

All simulations of water near model silica surfaces were 
performed using a modified version of the DL_POLY soft¬ 
ware package (version 2.18)P 3 I, following the work of Hu 
and Weekd^! Simulation cells 45.6 x 43.9 x 180 A 3 
in volume were used to simulate 2468 SPC/EP 5 water 
molecules at model silica surfaces with a buffering liquid- 
vapor interface far from the surface, to maintain a con¬ 
stant coexistence pressure while in the canonical ensem¬ 
ble (The large size of the cell in the ^-dimension is also 
useful for evaluating electrostatics in slab like geome¬ 
tries.). A temperature of 298 K was maintained using 
a Berendsen thermostat!^. Lennard-Jones interactions 
were truncated and shifted at a distance of 11 A. The 
slab corrected Ewald summation method of Yeh and 
Berkowitz 57 was used to handle the evaluation of elec¬ 
trostatic interactions, with a real space cutoff of 11 A, 
switching parameter a = 0.3 A -1 , and a maximum num¬ 
ber of k -space vectors of k x = 10, k y = 10, and k z = 30 
in the x —, y—, and z—directions, respectively, such that 
the k -space sums run from — k x to k x , for example. For 
the patterned dipolar surfaces discussed in Figure [6j a 


larger number of A;-vectors is needed, such that k x = 20, 
k y = 20, and k z = 30 for these simulations. 

B. Protein Simulation 

Simulations of HFBII (PDB ID: 2B97^P in SPC/EP 
water were performed using the GROMACS 4.5.3 soft¬ 
ware packagd^ and the AMBER 94 force fielcP^, follow¬ 
ing the work of Patel and Garde 18 . As with the silica 
surfaces, simulations were performed in the presence of 
a buffering liquid-vapor interface, now with one in each 
^-direction. The protein atoms were held fixed at the 
center of the liquid slab. Production simulations used 
for analysis totaled 6 ns in length, and a constant tem¬ 
perature of 300 K was maintained using the canonical 
velocity rescaling algorithm^. All Lennard-Jones inter¬ 
actions were truncated at a distance of 10 A, and long- 
ranged corrections to the energy and pressure were ig¬ 
nored. Electrostatic interactions were evaluated using 
the particle mesh Ewald method 61 . 

ACKNOWLEDGMENTS 

This work was supported by the National Science 
Foundation (Grants CHE0848574 and CHE1300993). 
We are grateful to Jocelyn Rodgers for helpful discus¬ 
sions. RCR also acknowledges Amish Patel for stimu¬ 
lating discussions and Erte Xi for assistance with the 
generation of the protein interface. 

Appendix A: Multipole Moment Expansion of 
Gaussian-Smoothed Charge Densities 

Here, we present the multipole expansion for the 
Gaussian-smoothed charge density, p qa ( r), and how it 
relates to that of the bare charge density p q { r). For gen¬ 
erality, we consider the interaction energy between two 
d-dimensional Gaussian-smoothed charge densities, with 
centers and r j , where the position vector is defined by 

r = (xi,X 2 , ...,x d ) , (Al) 

and d £ N. We consider the interaction energy of the 
two charge distributions under consideration, p qa { r^) and 
pY ( r j ), respectively, to be of the form 

w{Tij) = J dr J dr'pf(r - r- r 3 -) ^ ^ ^ , 

(A2) 

where we use the 3-dimensional Coulomb potential. How¬ 
ever, the multipoles themselves depend only on the di¬ 
mensionality of the charge density, and the relationship 
between the multipoles of the bare and smoothed charge 
densities described below are independent of the form of 
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the potential. This interaction energy can be rewritten In general, the nth order multipole moment A4 a (n) is 
as a k-space integral, given by 


w ( r v) = J ^r(-k)pf( k )e ik ' ry (A3) 

where we have defined the d-dimensional Fourier trans¬ 
form and inverse Fourier transform of a function / as 

/(k) = J dre~ ik ' r f(r) 

and 

/w =^/ dke *' /(k) - 

respectively. 

Now to examine the asymptotic behavior as k -> 0, 
we Tayor expand the smoothed charge densities about 
k = 0, such that 

pT( k) = £ T k n ; . v-pHO), (A4) 

. 1 ’ 


where Vk is the d-dimensional gradient with respect to 
k. We can then insert lA4l into lA3l to obtain 

IK'pTio) • Hv r H 

I to •I in • 

•[V”Vf( 0).(iV P )"*]i (A5) 

Now, we can define 

l ^VZpr(0)mMUni), (A6) 

Tli'. 

such that 

MU»i) =—, [ drpf(r)r”% (A7) 

np. J 


and A4 a (n) is the nth multipole moment of the smoothed 
charge distribution. Finally, we can rewrite A5 


as 


w(r«) = E \MKrn) • (-V P )"*] • [M?(nj) ■ V?'] \ 

n i 5 n 7 - ^ 

(AS) 

so that the energy is now expressed in terms of the mul- 
tipole moments of the smoothed charge distributions. 

One may then inquire into how these multipole mo¬ 
ments relate to those of the bare charge densities, p q ( r). 
In order to relate the two sets of multipoles, we first 
consider the Fourier transform of the smoothed charge 
density, which, using the convolution theorem, can be 
written as 


p qa (k) = p q (k)p G (k), (A9) 

where 

/3 G (k) = e -iL T“. (A10) 


m =0 


M*(n) = ^ E (I)^ (n " m) (°) ® 4 m) (°)» ( A11 ) 


where f^(0 ) = V£f( k) 


is a tensor of rank n, 


k=0 


indicates a symmetric outer product, and 


n\ 


m\(n — mn)\ 


is the binomial coefficient. The gradients of the k-space 
Gaussian function are given by 


4 n) (k) = (-l)V^/ 4 H n (y), 


kcr 


(A12) 


such that H n (ax) is a rank n tensor-analog of the Her- 
mite functions with elements 


Hij... v (ax;n) = (-l) n e a 


d n 


dxidxn • • • dx , 


(e- a2x2 ) , 


(A13) 

where a is a constant and x = (aq, aq, •••, xj) is a general 
d-dimensional vector. 

All odd derivatives of pc( k) will vanish at k = 0 due 
to symmetry, therefore, we can rewrite All 


as 


M*(n) = l - E (-l) m (j\P q{n ~ m \ 0) ® A - (|) 


m =0 
raEE 

(A14) 

where A m (a) = [H m (ak)] | fc _ Q and E is the set of even 
whole numbers. 

Equation |A14 can be written in the equivalent form 

Am ( 2 )’ 


71 / -« \ yyi • 

AA a (n) = M(n) + ^ - — 1 M(n — m) 


m =2 
raEE 


m! 


(A15) 

making the relation between A4 a (n) and A4(n ) apparent. 
Therefore, in order for A4 a (n) = A4(n) to hold, where 
A4(n) is the nth multipole moment of the bare charge 
distribution p q , all multipoles of the bare charge density 
p q of order less than n — 1 and of even (odd) order, for n 
even (odd), must be identically zero, 


M a (n) =M(n ) 
where 


M(s) = 0Vs = n-Z, (A16) 


/ 2,4,6, for n even , . 

| 2,4,6, ...,n — 1; for n odd. ' ' 

To illustrate this condition, we present the first few 
multipole moments of the Gaussian-smoothed charge 
density. The monopole moment of p qcr is trivially given 
by A4 a (0) = A4(0), and note that for neutral charge 
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distributions the monopole moment is zero. The dipole 
moment, n = 1, is also trivially given by 

M a (l) = M(1), 

illustrating that the dipole moment of p q is conserved 
upon Gaussian-smoothing. In addition, the quadrupole 
moment is given by 

AT(2) = M(2 ) - l -M{ 0) ® A 2 (|) 

2 

= M(2) + ^(0)I 3 , (A18) 

while the Gaussian smoothed octupole is similarly given 
by 

2 

M a (3) = 7W(3) + ^-Vf(l)®I 3 , (A19) 

where I n is the n x n identity matrix. For neutral 
charge distributions, like non-ionic molecular charge dis¬ 
tributions (water), both the dipole and quadrupole mo¬ 
ments are conserved upon Gaussian-smoothing, but not 
higher order moments. From the work of Wilson and 
PrattPSl, it follows that only the potential difference 
across two phases in a slab-like geometry is conserved 
upon Gaussian-smoothing of the charge density (or the 
potential), because this depends only on the dipole and 
quadrupole moments of a molecular fluid. However, con¬ 
trary to a previous report^, Vs (z) ^ V(^), because 
higher order multipoles are modified by smoothing. 
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